Corregionlization
Semipositive definite matrix spherical model
mat1 <- cbind(c(30, 30, 30),
c(30, 50, 30),
c(30, 30, 35))
#matriz definida positiva "cercana"
mat1 <- data.frame(as.matrix(nearPD(mat1)$mat))
names(mat1) <- c("NO2", "O3", "NOX")
row.names(mat1) <- c("NO2", "O3", "NOX")
pander::pander(mat1)
| NO2 |
30 |
30 |
30 |
| O3 |
30 |
50 |
30 |
| NOX |
30 |
30 |
35 |
Semipositive definite matrix spherical hol.
mat2 <- cbind(c(13.02, 24.5, 18.739),
c(24.58, 46.4, 35.36),
c(18.73, 35.36, 26.95))
mat2 <- data.frame(as.matrix(nearPD(mat2)$mat))
names(mat2) <- c("NO2", "O3", "NOX")
row.names(mat2) <- c("NO2", "O3", "NOX")
pander::pander(mat2)
| NO2 |
13.02 |
24.54 |
18.73 |
| O3 |
24.54 |
46.4 |
35.36 |
| NOX |
18.73 |
35.36 |
26.96 |
G_stat object
Semivariogram
vgmno2 <- vgm(psill = mat1[1, 1],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 1],
model = "Hol",
range = 2294))
vgmo3 <- vgm(psill = mat1[2, 2],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 2],
model = "Hol",
range = 2294))
vgmnox <- vgm(psill = mat1[3, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[3, 3],
model = "Hol",
range = 2294))
Cross semivarogram
vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 2],
model = "Hol",
range = 2294))
vgmno2_nox <- vgm(psill = mat1[1, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 3],
model = "Hol",
range = 2294))
vgmno3_nox <- vgm(psill = mat1[2, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 3],
model = "Hol",
range = 2294))
Gstat
remove_na <- function(frame, vari_) {
# Remove na from sp object
datos1 <- frame
bool <- !is.na(datos1@data[vari_])
datos1@data <- datos1@data[bool, ]
datos1@coords <- datos1@coords[bool, ]
return(datos1)
}
coordinates(datos) <- ~ X + Y
g_st <- gstat(NULL,
id = "NO2",
formula = NO2 ~ X + Y,
model = vgmno2,
data = remove_na(datos, "NO2"))
g_st <- gstat(g_st,
id = "O3",
formula = O3 ~ Y,
model = vgmo3,
data = remove_na(datos, "O3"))
g_st <- gstat(g_st,
id = "NOX",
formula = NOX ~ Y,
model = vgmnox,
data = remove_na(datos, "NOX"))
#cross
g_st <- gstat(g_st,
id = c("NO2", "O3"),
model = vgmno2_o3)
g_st <- gstat(g_st,
id = c("NO2", "NOX"),
model = vgmno2_nox)
g_st <- gstat(g_st,
id = c("O3", "NOX"),
model = vgmno3_nox)
pander::pander(do.call(rbind, g_st$model)[, 1:3])
| NO2.1 |
Hol |
13.02 |
2294 |
| NO2.2 |
Sph |
30 |
6096 |
| O3.1 |
Hol |
46.4 |
2294 |
| O3.2 |
Sph |
50 |
6096 |
| NOX.1 |
Hol |
26.96 |
2294 |
| NOX.2 |
Sph |
35 |
6096 |
| NO2.O3.1 |
Hol |
24.54 |
2294 |
| NO2.O3.2 |
Sph |
30 |
6096 |
| NO2.NOX.1 |
Hol |
18.73 |
2294 |
| NO2.NOX.2 |
Sph |
30 |
6096 |
| O3.NOX.1 |
Hol |
35.36 |
2294 |
| O3.NOX.2 |
Sph |
30 |
6096 |
Plot model model and estimation
plot(variogram(g_st),
model = g_st$model,
pl = T)

Prediction plot.
prediction_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 100000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
geom_tile() +
scale_fill_viridis_c() +
theme_void()
return(plot)
}
variance_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
var <- paste(variable, ".var", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
geom_tile() +
scale_fill_viridis_c(option = "inferno",
direction = -1) +
theme_void()
return(plot)
}
cv_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
var <- paste(variable, ".var", sep = "")
aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
geom_tile() +
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
return(plot)
}
pl1 <- prediction_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl2 <- variance_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl3 <- cv_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/juan/Documents/Inves/Webinar/Day1/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
ggplotly(pl1)